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ABSTRACT 


Maximum range and minimum time maneuvers for a roeket are vital to the war 
fighter since weapons that can intercept targets faster or travel farther provide advantages 
that increase the chances of mission success. Yet, guidance laws for tactical rockets are 
based on principles that have been in use since World War II. This thesis applies 
Pseudospectral optimal control theory to assess the achievable performance of a guided 
rocket in the terminal stage of flight. The performance is compared to a baseline obtained 
using conventional control (proportional navigation). For a fictitious missile, an optimal 
control solution is shown to increase the range of the missile by nearly 300% in the mid- 
to-terminal phase of flight. A minimum time solution showed a reduction in flight time 
by about 40% compared to the conventional proportional navigation-based approach. 
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I. 


INTRODUCTION 


A. BACKGROUND 

One area that needs to be addressed prior to any diseussion is: what is the 
difference between a missile and a rocket? While the exact definition varies from source 
to source, the general consensus throughout the aerospace community is that “a missile 
has a guidance system or brain to get it to its destination and a rocket just goes where it 
is initially pointed” [1]. The U.S. Bureau of Naval Personnel defines a guided missile as 
“an unmanned vehicle that travels above the earth’s surface; it carries an explosive war 
head or other useful payload; and it contains within itself some means for controlling 
its own trajectory or flight path” [2], The terms guided rocket and missile are used 
interchangeably throughout this thesis. 

Prior to World War I, missiles and bombs were essentially shot or dropped on a 
certain trajectory, and there was no capability to control the device after its release. The 
concept for guided missiles was first conceived during World War I. The tactical use of 
aircraft in World War I raised the question as to whether a plane could be remotely 
controlled and used to bomb a target. The first guided missiles that emerged were small 
versions of the aircraft that were currently being flown by utilizing radio control [3]. It 
was not until World War II though, that tactical guided missiles began to emerge on the 
scene. 

Today, guided missiles are an integrated part of the United States Order of Battle. 
The guidance system for a missile allows it to find a path to intercept the target “in spite 
of initial launching errors, in spite of wind or other forces acting on the missile, and in 
spite of any evasive actions the target may take” [2]. There are currently several types of 
guidance modes that are used in guided missiles, including command guidance, homing 
guidance, and navigational guidance. A missile may have more than one type of guidance 
system, with different systems engaging at different stages of flight. 
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Early Guided Rockets 


The first guided missiles appeared in Germany during World War II. The most 
widely known missiles that were produeed during this time were the V-1 and the V-2. 
The V-1 was a “pilotless, pulsejet, midwing, monoplane, laeking ailerons but using 
eonventional airframe and tail eonstruetion” [2], shown in Figure 1. The guidanee and 
eontrol for the missile was done using an internally stabilized gyro and the existing 
eompass guidanee. These “missile fixed rate gyros were used to provide attitude rate 
signals for stabilization” [4]. The missile was not very aeeurate but its massive warhead 
eaused significant damage if it hit the intended target. However, the V-1 was slow and 
approximately 95% of the missiles launched against England were brought down by anti¬ 
aircraft fire [2]. 



Figure 1. V-1 Missile 

Source [5]: Warbirds Resource Group, “Fi-103/V-l ‘Buzz Bomb’,” 2015 (November 12). 

The V-2, shown in Figure 2, was a much larger than the V-1 and had a liquid-fuel 
rocket propulsion system. The missile was launched using a vertical launch system. Prior 
to launch the system was loaded with a command angle that would cause the missile to 
rotate to a predetermined angle, between 41- to 47-degrees, shortly after launch. The 
launch speed required to reach the target would also be determined prior to launch and 
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when the missile reaehed the desired veloeity, a radio eommand would be sent to the 
roeket and the propulsion system would be shut down. This was aeeomplished by 
utilizing the Doppler Effeet to determine the missile veloeity. The missile utilized a Wien 
bridge and the beat frequeney of the return signal would be sent to the bridge in addition 
to the ground station’s frequeney seeond harmonie. Prior to launeh a beat frequeney 
would have been ealeulated eorresponding to the desired eutoff veloeity. When the 
reeeived frequeney matehed the preset frequency, a command was transmitted to the 
missile to initiate the cutoff of the propulsion system [4]. Once the propulsion system was 
shut down, the range of the rocket was only about 320 km, it would reach speeds of 
approximately 1475 m/s. While the missile was not very accurate, its high speed made 
countermeasures impossible and “no V-2 missile was ever intercepted or shot down by 
anti-aircraft fire” [2]. 



Figure 2. German V-2 Rocket 

Source [6]: Arkell, H., 2014, “Death from Above without Warning: 70 Years After the 
First One Fell, Interactive Map Reveals just Where Hitler’s V2 Rockets Killed 
Thousands of British Civilians in Final Months of WW2,” 2015 (November 12). Daily 
Mail. 


On the other side of the Atlantic, the United States initiated the development of 


the Lark missile, shown in Figure 3, under a guided missile development program in 
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response to the effeetive Kamikaze attack against U.S. ships. In December of 1950, the 
first successful intercept of a guided missile against an unmanned aircraft was made [7]. 
The missile used an active continuous wave radar onboard the missile. Reflections from 
the target would be displaced in frequency by the Doppler Effect. However, the missile 
would also be receiving inputs from the ground, and once the stationary inputs were 
rejected, the missile would be left with the corresponding frequency to the moving 
target [8]. However, this type of guidance presented several fundamental problems. The 
first problem was that the Doppler shift from the target was not very large. The second 
problem was that there was a large amount of clutter in the return signal. The need for 
low-noise microwave sources became a critical part of the continuous wave radar 
development as a result of these issues [8]. At the end of World War II, the inertial 
sensors were limited in their accuracy due to the fiction in the ball and pivot bearings 
causing errors. Later, by utilizing cylindrical gas bearings that were externally 
pressurized, the accuracy of the internal sensors increased significantly [4]. These 
improvements paved the way for modern guidance systems. 



Figure 3. Lark Guided Missile 

Source [9]: Parsch, A., “Lark,” Encyclopedia Astronautica, 2015(November 14). 
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2. Modern Missile Guidance, Navigation, and Control Systems 

A basic block diagram for a modern guidance, navigation, and control (GNC) 
system for a missile is shown in Figure 4. The main components include a weapon 
control system, sensors, flight controller, and airframe subsystems that control the 
missile. A more in-depth description of the GNC system is found in [10]. The weapon 
control system is ground or air based and allows the user to send commands to the missile 
and also receives flight data that is sent from the missile flight computer. 


Worid data 



Figure 4. Basic GNC System for a Missile 

Source [10]: Hawley, P. A., and Blauwkamp, R. A., 2010, “Six-Degree-of-Freedom 
Digital Simulations for Missile Guidance, Navigation, and Control,” John Hopkins APL 
Technical Digest, Vol 29(No. 1) pp. 71-84. 


The missile also contains two types of sensors, inertial sensors for determining 
where the missile is located in space and terminal sensors for determining where the 
target is in relation to the missile. The most commonly used inertial sensors include 
inertial measurement units, inertial reference units, and GPS. One or more of these 
sensors may be used to measure and/or help estimate the states of the missile. The 
terminal sensors for a missile will be based on the type of guidance system the missile 
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utilizes. Some of the more common sensors include radar, radio frequency seekers, and 
infrared seekers. 

The flight controller is the brains of the missile and includes the navigation 
processor, guidance processor, and autopilot. The navigation processor uses external 
commands from the weapons control system and inertial sensor data to estimate the states 
and rates for the missile. This information is passed on to the guidance processor. The 
guidance processor receives inputs from the terminal sensor, the navigation processor and 
the weapons control system. It then computes and outputs a guidance command that is 
sent to the autopilot for implementation and sends sensor commands to the terminal 
sensors. The guidance law for a missile is implemented in the guidance processor, while 
the autopilot is used to convert the guidance signal into flight commands that will be sent 
to the airframe [10]. 

The flight commands from the autopilot are sent to either the propulsion system 
or the flight control actuators for the missile. The propulsion system will make any 
needed adjustments and the system may also output the mass of the propulsion that has 
been burned to determine the airframe mass changes. The flight control actuators adjust 
the flight control surfaces to produce an aerodynamic response. The propulsion forces 
and moments, airframe mass changes, and aerodynamic forces and moments are used to 
estimate the state and state rates for the missile. These states and state rates are used 
along with inertial sensors and environmental data to determine where the missile is 
located in space. 

The U.S. Armed Forces currently has dozens of different types of guided 
rockets in their Order of Battle today that can be launched from land, sea, and the air. 
They range from ballistic guided missiles (e.g.. Trident II D5), shown in Figure 5, 
to laser-guided rockets (e.g., APKWSTM), shown in Figure 6, to active radar homing 
(e.g. AIM-120 AMRAAM), anti-radiation homing (e.g., AGM-88 HARM), and GPS 
guided (e.g., UGM-133 Trident II). While tactical missiles that are used today are able to 
hit targets much more accurately, this is mainly due to improvements of the sensors in the 
missiles rather than the evolution of the overarching principles governing the guidance 
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control. In fact, the same guidance principals that were used for the Lark in 1950, are still 
being used today [7]. 



Figure 5. Trident II D5 Fleet Ballistic Missile 

Source [11]: Lockheed Martin, 2014, “Trident II D5 Missile Reaches 150 Successful Test 
Flights,” 2015(Nov 13). 



Figure 6. Laser-Guided APKWS Rocket 

Source [12]: Sax, L. R., 2015, “U.S. Army Acquires APKWS™ Laser-Guided Rockets 
for Immediate Deployment,” 2015(Nov 13). 


For example, consider the AGM-114A is the original Hellfire missile that was 


used by the Army. It was fired from the air and one of the delivery modes for the missile 
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was Lock-on After Launch (LOAL). The LOAL profile allowed the airerew to fire a 
missile in the general direction of a target without having a line-of-sight look on the 
target. Within this profile, there are two delivery modes where the airoraft is able to be 
masked behind terrain when firing to reduoe the risk to the aircraft: low and high. In 
LOAL-LO, shown in Figure 7, the missile is to olear 79.25 m (260 ft) of terrain with the 
aircraft located 600 m behind the highest terrain point. In LOAL-HI, shown in Figure 8, 
the missile oan olear 304.8 m (1000 ft) of terrain with a 1500 m airoraft standoff distanoe. 
The missile utilizes a laser seeker looated in the missile nose, whioh is programed prior to 
launch, to identify a speoifio laser oode and look onto the laser emission. For this missile, 
the laser seeker is the terminal sensor and after the missile obtains a look on the target, 
the information from the seeker is sent to the guidanoe prooessor whioh then determines 
the appropriate steering command to be sent to the flight oontrols and the missile then 
tracks to the target [13]. While the laser seeker allows for inoreased aocuracy in hitting 
the target, oompared to the missiles used in World War II, the guidanoe law for the 
missile remains unohanged. 



LOAL-LO Trajectories (AGM 114-A) 


Figure 7. AGM-114A LOAL-LO Trajeotory 

Source [13]: Lange, A. W., 1998, "Hellfire: Getting the most from a Lethal Missile 
System," Armor, Jan-Feb, pp. 25-30. 
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LOAL-HI Trajectories (AGM 114-A) 

Figure 8. AGM 114-A LOAL-HI Trajectory 

Source [13]: Lange, A. W., 1998, "Hellfire: Getting the most from a Lethal Missile 
System," Armor, Jan-Feb, pp. 25-30. 


B. MOTIVATION AND SCOPE 

A paper written by Lin [14], states: 

Many long- and medium-range missile guidance studies have shown that 
optimal trajectory shaping promises an extended range with more 
favorable end-game conditions... However, in a three-dimensional, target- 
intercept flight, direct application of the optimal control theory will result 
in a two-point boundary value problem that involves several arbitrary 
parameters so that analytical solutions cannot be obtained without a lot of 
approximations. The problem can be further complicated with the lift, 
thrust and drag, and control constraints forced by structural and angle-of- 
attack limits. This increases the computational time so that it is not 
feasible to implement the resulting solutions of the missile performance 
within an onboard digital computer. In view of this complexity in the 
problem setup, either indirect methods or direct methods based on 
nonlinear programming are used to solve the sensitivity and convergence 
problems. However, both methods require very fast onboard 
microprocessor technology for real-time, on-line operation. [14] 
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It is now almost 30 years since this article was written, and the same challenges 
still exist today. The dynamic equations that govern a three-dimensional 6DOF model for 
a missile are complex. The equations consist of nonlinear trigonometric functions that 
present potential singularity problems. In addition, aerodynamic coefficients typically 
come in tables aerodynamic forces must be “built-up” so that the dynamics can be 
evaluated. 

The boundary value problem for solving the optimal control contains a set of 
prescribed boundary conditions and also a set of constraints on the problem variables. 
The boundary conditions can be set values or functions. Airframes have limits that are 
applied to the angle of attack or acceleration, for example. The solution for the boundary 
value problem must therefore be constrained by these limits. In order to solve a boundary 
value problem, the underlying differential equations (e.g., the 6-DOF dynamics) need to 
be solved while simultaneously satisfying the constraints of the bounds [15]. The large 
number of constraints that are applied to the guidance of a missile, as described in the 
quote above, have led to difficulties in solving these boundary value problems in the past. 
Today, however, advanced tools are available to overcome these challenges, [15], [16]. 

As mentioned above, the principles that are used to control guided rockets today 
are the same principles used on some of the first guided rockets. The motivation behind 
this thesis is to examine how advances that have been made in optimal control theory 
might be used to change the way that guided rockets are controlled. The goal is to 
determine if a tactical advantage, such as enhanced range, can be obtained via a software 
solution only. Such a solution is advantageous in that its implementation would not 
change any of the physical characteristics of the missile. 

In a warfighting environment, the difference between life and death for a solider 
comes down to being able to maintain a tactical advantage over the enemy. This 
advantage can come in many forms. For example, if applying an optimal control solution 
to the guidance of a rocket can increase the range or reduce the time of flight, it can 
provide our troops with an increased advantage in the battlefield. 
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The basic principles of the GNC that are used on tactical missiles, not only apply 
to missiles within the Earth’s atmosphere but also apply to other systems such as 
intercontinental ballistic missiles, spacecraft performing close proximity rendezvous 
operations, and UAVs. Therefore, the results of this thesis could influence how guidance 
is done for these systems as well. 

C. THESIS ORGANIZATION 

The remainder of this thesis is organized as follows. Chapter II introduces the 
model of the fictitious rocket used in this work. In addition, the 3DOF and 6DOF models 
that can be used to simulate the flight of the guided rocket are introduced and the pros 
and cons of each are discussed. The software used to determine the coefficients of lift and 
drag is introduced, as well as the development of the coefficient of lift and drag 
equations. Chapter III presents the conventional proportional navigation approach that 
has long been used for guidance of the rockets. In addition, a maximum range maneuver 
is simulated to establish a baseline performance that can be used to compare against the 
results achieved through optimal control. In Chapter IV the Optimal Control Theory is 
introduced in order to provide the necessary background on the subject for the reader. In 
addition, this chapter introduces the optimal control software DIDO. This software 
package is used to solve the optimal missile trajectories. Chapter V discusses the 
maximum range problem and the steps that were taken in order to solve it. In addition, 
some of the challenges that were encountered when attempting to obtain the solution are 
discussed. Chapter VI presents some conclusions and suggestions for future work. 
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II. GUIDED ROCKET MODEL 


This chapter covers the development of the model for the tactical missile that will 
be used for simulations in this thesis. It introduces not only the physical characteristics of 
the model but also the assumptions that were made in regard to the model. The reference 
frames and equations of motion for the missile are also introduced. This chapter also 
shows the development of the coefficients of lift and drag for the missile using the 
Missile DATCOM software [17]. 

A. DEVELOPING THE MODEL 

The first step in creating the model was to choose the physical configuration of a 
fictitious missile. In order to create a process that could be adapted to any missile, a 
generic rocket was created that is similar to the larger tactical missiles currently fired 
from aircraft, which are categorized as medium-range missiles. This size of missile was 
chosen with the intent that the model could later be scaled up or down in order to mimic 
different systems. 

1. Guided Rocket Configuration 

a. Physical Characteristics 

The model for the rocket, shown in Figure 9, has a cylindrical body with a cone 
shaped nose and contains four wings and four tail fins. 
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Figure 9. Image of a Fietitious Guided Missile 


The fictitious guided rocket is 4.14 m in length with a wingspan of 1.01 m. The 
diameter of the rocket is 0.254 m around the tube, shown in Figure 10. The center of 
gravity for the guided rocket was placed at 2.08 m. The mass of the guided rocket 
is 360 kg. The guided rocket has two different sets of fins attached. The first set is a set of 
tail fins that are located on the tail end of the rocket and the second set is in the middle of 
the rocket and are referred to as the wings. 



Figure 10. Image of the Model from the Side with Dimensions 
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Each wing of the guided roeket is 0.432 m with a width of 0.216 m where the 
wing attaehes to the rocket body. The leading edge of the wing is tapered and has 
deereasing slope for the first two thirds of the wing span and then levels out for the last 
third. The width of the wing at the tail end is 0.108 m as shown in Figure 11. There are 
four wings on the missile that are mounted at 45° increments around the missile body 
(see Figure 12). 



Figure 11. Missile Wing Dimensions 
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In addition, there are four tail fins that are mounted on the tail end of the missile. 
The fins are 0.178 m with a width of 0.216 m where the fin attaehes to the body. The fins 
are tapered and at the outside edge are a width of 0.144 m, shown in Figure 13. 



Figure 13. Missile Fins 

The 6DOF equations of motion whieh are discussed in the next section of this 
chapter, also require a reference aerodynamic chord and reference aerodynamic span. 
These values are 0.063 m and 1.12 m, respectively. A summary of the physical 
characteristics of the missile is given in Table 1. 
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Table 1. Summary of Physical Characteristics for a Fictitious Missile 


Parameter 

Symbol 

Value 

Units 

Mass 

m 

360 

kg 

Length 

1 

4.14 

m 

Diameter 

d 

0.83 

m 

Wing span 

Iws 

1.01 

m 

Fin span 

^fs 

0.61 

m 

Center of gravity 

CG 

2.08 

m 

Referenee aerodynamic chord 

c 

0.06 

m 

Reference aerodynamic span 

b 

1.12 

m 

Cross sectional area 

^ref 

0.05 

m' 


b. Performance Characteristics 

In addition to designing the missile geometry, performance characteristics need to 
be created for the missile in order to be able to carry out control system design. These 
numbers were estimated by looking at a variety existing platforms and estimating values 
that would be in the same order of magnitude for a similar model. When the process 
described in this thesis is applied to an actual missile and more accurate numbers are 
available, the results that are produced from the control design and analysis will better 
reflect reality. The estimated performance values are listed in Table 2. 


Table 2. Performance Characteristics for a Fictitious Missile 


Parameter 

Symbol 

Value 

Units 

Minimum velocity 

Vmin 

100 

m/s 

Maximum velocity 

Vfnax 

639 

m/s 

Minimum thrust 

T 

^ min 

0.00 

N 

Maximum thrust 

T 

^ max 

10000 

N 
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2. 6DOF Dynamics Model 

In a 6DOF dynamic model, all six degrees of freedom for the guided rocket are 
taken into account in regard to the rocket motion. There are three translational degrees of 
freedom and three rotational. The six degrees of freedom result in a total of twelve non¬ 
linear state equations. 

a. Reference Frame 

The 6DOF model uses a flat-Earth NED as a reference frame for position, as 
shown in Figure 14. 



Figure 14. Body-Axis Reference Frame for Missile 

All other states use a body-reference frame, as shown in Figures 15 and 16. The 
flight path angle, y, illustrated in Figure 15, is the angle between the local horizon and the 
velocity vector. The angle of attack, a, is the angle between the velocity vector and the x- 
axis in the body-reference frame. The pitch angle, 6, is the angle between the local 
horizon and the x-axis in the body frame. 
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Figure 15. Flight Angles in the X-Z Body Reference Frame 

The side-slip angle, fi, shown in Figure 16, is the angle between the velocity 
vector and the x-axis in the body-reference frame. Yaw is the rotation of the x-y body- 
reference frame about the z-axis. The yaw angle, )//, is the angle between the x-axis in the 
body-reference frame from a locally referenced position at the beginning of the rotation. 



y 

Local Reference 


Figure 16. Flight Angles in the X-Y Body Reference Frame 
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The heading angle, a, is the angle between the velocity vector and the x-axis in 
the x-y plane, where the velocity vector is in the direction of travel from the COG 
through the nose of the missile, shown in Figure 17. 



Figure 17. Definition of Flight Path Heading Angle 

b. Dynamics Equations 

For a 6DOF model, the equations of motion are a function of system states, 
forces, moments, and a number of parameters. The variables that represent the moment, 
forces, velocities, and parameters are listed in Table 3, while a representative set of states 
for the 6DOF model are listed in Table 4. 
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Table 3. 6D0F Forces, Moments, and Velocity Components 



Roll (x-axis) 

Pitch (y-axis) 

Yaw (z-axis) 

Angular Rates 

P 

q 

r 

Velocity Components 

u 

V 

w 

Aerodynamic Force Components 

X 

Y 

Z 

Aerodynamic Moment Components 

L 

M 

N 

MOI about each axis 

^XX 

lyy 

^ZZ 


Table 4. 6DOF States for Guided Rocket 


Name 

Symbol 

Units 

Velocity 

V 

m/s 

Angle of Attack 

a 

rad 

Sideslip Angle 

P 

rad 

Roll rate 

P 

rad/sec 

Pitch rate 

q 

rad/sec 

Yaw rate 

r 

rad/sec 

Roll angle 

¥ 

rad 

Pitch angle 

e 

rad 

Yaw angle 


rad 

X-coordinate (Earth Axes) 

xe 

m 

Y-coordinate (Earth Axes) 

ye 

m 

Altitude 

ze 

m 


For each of the states that are listed in Table 4, there is an associated dynamic 
equation. The standard 6DOF equations of motion from [18], [19] include the position 
equations (2.1), the control rate equations (2.2), the angular velocity equations (2.3), the 
Euler angle rate equations (2.4), and the force and moment equations (2.5). 


21 


























= (cos 6 cos Vcos a cos P + (sin (j) sin 6 cos y/ - cos (j) sin y/) Vsin p 
+(sin sin ^ + cos (j) sin 6 cos y/) V sin a cos (3 

= (cos 6 sin y/) V cos a cos f3 + (cos (j) cos ^ + sin sin 6 sin Vsin (2.1) 

+(cos (j) sin 6* sin ^ - sin (j) cos Vsin a cos P 

= (- sin Vcos a cos P + (sin (j) cos 6) Vsin p + (cos (/) sin 6*) Vsin a cos P 


X Y 

V = (- g sin 6) cos a cos P + {— + g sin ^zi cos 6) sin P 
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+(— + g cos (j) cos 6) sin a cos P 
m 

1 g 

d = -(Z cos « - X sin «) H-^— (cos S cos 6 cos « + sin ^ sin a) 

mVco?,p VcosP 
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mV 
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'(p (^ZZ ^YY^V^XZ^P) 
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^ XX ^zz ^xz 


-(N+(/^ IjY^^P ^XzV^ 


— (M-(Ijo,-Izz)pr-Ixz(P -r )) 


^xx^zz ^ xz 


- {L + (/yy 7^2 )qr + I^qp) 
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6 = qcos(/)-rsm(/) 
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(2.3) 


(2.4) 
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X = i p(z)F%^^C^ {a,p,M,S) + T^ 

Y = ^p{z)V^S^^^C,{P,MA) + T^ 

Z = 1 p(z)FX^^C, {a,p,M,5) + T^ 

L = ^p{z)F^S^^^bC,{P,M,5) 

M = ^p{z)F^S^^^cC^{a,p,M,5) 

N = ^p{z)F^S,^^bCSP,MJ,5) 


(2.5) 


While a 6D0F model more aeeurately represents the performanee of a missile, a 
3DOF model is often suffieient for optimal trajectory design [18]. 

3. 3DOF Dynamics Model 

In order to simplify analysis, a 3DOF model is often used in lieu of the 6DOF 
equations of motion. The 6DOF model involves multiple assumptions. These 
assumptions necessarily hold true in the 3DOF model as well. Additional assumptions are 
made that set the necessary conditions where the 3DOF equations of motion represent the 
6DOF model. All assumptions for the 3DOF model are described in the next section. 

a. Model Assumptions 

The first assumption that is made (in both the 6DOF and 3 DOF models) is that the 
Earth is flat and is non-rotating. By assuming that the Earth is non-rotating, the Coriolis 
Effect that is caused by the rotating earth goes to zero. This model also makes the 
assumption that the missile is not traveling far enough of a distance for round-earth 
effects to be taken into account. Due to the fact that the Earth is not a perfect sphere, 
there are variations in the gravitational constant based on where you are located on the 
surface of the Earth. The model assumes that there is no variations in the gravity so the 
standard constant value is used. 
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It is also assumed that the missile has a eonstant mass. In the flight of an aetual 
missile, the mass will ehange through the beginning stages of flight as the fuel is burned. 
For the problem being examined in this thesis, the missile is in the mid-to-terminal phase 
of flight and it is assumed that all the propellant has been burned. When the entire missile 
trajeetory is taken into aeeount, mass variations need to be eonsidered. The missile is also 
treated as a rigid body vehiele. By assuming that the missile is a rigid body rather than a 
flexible body, it eliminates the distortional effeets that oeeur due to the additional forees 
and moments eause by aeroelastieity in the strueture (e.g., flutter of the wings and fins). 
Additionally, the eross-produets of inertia are eonsidered to be negligible. 

Missiles are launehed in many environments and regions in the world. For most 
areas over land, the ground is loeated above the sea level and when the missile impaets 
the ground is it not at a mean sea level of zero. The 6DOF and 3DOF models both 
assume that when z = 0, the missile is at sea level, in a standard atmosphere. In addition 
sinee the missile is symmetrie about the roll axis, the roll rate is assumed to be fixed at 
zero. Side-slip results from side foree that is produeed from an unequal lift on the 
airframe wings from either untrimmed flight surfaees or relative wind [20]. For this 
model is also assumed that there is no sideslip. The assumptions that were made for the 
model are summarized in Table 5. 
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Tables. Model Assumptions 
Flat, non-rotating Earth 
No wind 

No gravity variations 
No Coriolis effect 
Rigid body vehicle 
Constant mass 

Negligible cross-products of inertia 
Z = 0 is at sea level, standard atmosphere 
No side-slip 

Missile is symmetric about the roll axis 


b. Dynamic Equations 

For the 3DOF model dynamic equations, the variables that are used are defined in 
Table 6. 


Table 6. 3DOF Dynamic Equation Variables 


Name 

Symbol 

Units 

Down-range distance 

X 

m 

Cross-range distance 

y 

m 

Altitude 

z 

m 

Velocity 

V 

m/sec 

Heading angle 

a 

rad 

Flight path angle 

Y 

rad 

Angle of attack 

a 

rad 

Bank angle 

E 

rad 

Thrust 

T 

N 

Drag 

D 

N 

Lift 

L 

N 

Mass 

m 

kg 
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When the assumptions deseribed above are applied to the 6DOF equations of 
motion, the equations are redueed as follows: 


i = veosf cos cr 
y = vcosf sin cr 
z = vsinf 

v = —(T cos a-D)-g sin/ 
m 


& = —{T sina + Z) 
mv 


sin ju 
cos/ 


/ = — {T sin a + L) cos ju 
mv 


— cos/ 

V 


( 2 . 6 ) 


If it is assumed (as above) that the rocket motor is burned, thrust drops from the 
equations and Equation (2.6) can be further simplified. In this case the missile can be 
assumed to be in the mid-to-terminal phase of flight and the equations for the 3DOF 
model without the thrust terms reduce to: 


i = vcosf cos cr 
y = vcos/sincr 
z = vsm/ 

. -D 

V = - gsm/ 

m 

L sin jj, 

<j =- 

mv cos/ 

L g 

/ = —cos//- cos/ 

mv V 


(2.7) 


Equation (2.7) is the model that will be used in the remainder of this work. 


B. AERODYNAMIC MODEL 

1. Missile DATCOM Software 

The actual coefficients of lift and drag for a missile are determined through flight 
testing of the vehicle. However, since the values are needed for the design and simulation 
prior to actual testing, there is a need for software that predicts the aerodynamic 
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characteristics of a missile based on the design eonfiguration. There are several different 
models that ean be used in order to determine the eoeffieients of lift and drag for a 
missile. The software that was used for this thesis was the Missile DATCOM 
software [17]. The Missile DATCOM software is an aerodynamie design tool that allows 
users to input the physieal dimensions of a missile for a given eonfiguration and based on 
this, provides output of the required aerodynamie eoeffieients. 

The Missile DATCOM software is written in FORTRAN and allows the user to 
input the eharaeteristies of the missile that they are interested in, in the form of a 
eonfiguration file. The user is able to input the number of angles of attaek that that they 
desire and the Maeh number they are expeeting the missile to fly at. For the simplieity of 
this thesis, only one Maeh number was used. For a more eomplex and realistie model, a 
table of eoeffieients would need to be ereated in order to more aeeurately represent the 
aerodynamies. Dimensions of the missile are also needed whieh inelude the length of the 
nose seetion, diameter of the nose, number of fins and their loeation on the body, as well 
as dimensions of the fins. The program also allows for a trim setting to be applied and 
will output the assoeiated eoeffieients if the missile is able to be trimmed for the given 
angle of attaek. 

Onee the appropriate eonfiguration data is entered into the Missile DATCOM 
software, the program outputs a table of eoeffieients for a given set of trim settings and 
angles of bank. In the 6DOF model, the foree and moment equations. Equation (2.5), 
eaeh require an aerodynamie eoeffieient. There are six different eoeffieients, whieh are 
listed in Table 6. Of note, the roll moment eoeffieient, C/, should not be eonfused with 
the eoeffieient of lift. Cl. 
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Table 7. 6D0F Aerodynamic Coefficients 


Name 

Symbol 

Axial Force Coefficient 

Cx 

Side Force Coefficient 

Cy 

Normal Force Coefficient 

Cz 

Roll Moment Coefficient 

Cl 

Pitch Moment Coefficient 

Cn, 

Yaw Moment Coefficient 

c„ 


An example of one of the outputted aerodynamic coefficients is shown in Figure 
18. For ease of comparison all of the plots for the coefficient data have been normalized. 
Figure 18 shows the relationship between the angle of attack and the axial force 
coefficient. These numbers were all calculated at the assumed maximum velocity 
of 639 m/s (1.863 Mach). 
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Figure 18. Normalized Axial Force Coefficient versus Angle of Attack 
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When using the simplified 3D0F model, Equation (2.7), the equations of motion 
only require the values for lift and drag. The equations for lift and drag depend on the lift 
and drag eoeffieients. Cl and Cd as follows: 

i = t/OvX/Q 

2 ( 2 . 8 ) 

The eoeffieients of lift and drag ean be geometrieally determined from the axial 
and normal foree eoeffieients [18]. The resulting equations are: 


= -C^ cos a cos /^ - Cy sin/^- sin a cos (5 

Cy = Cjy sin a-Cy eos a 


(2.9) 


Sinee the 3DOF model assumes that there is no side-slip and that all turns are 
done in a eoordinated fashion where = 0, the equations further reduee to: 


Cjj - -Cyy eos a-Cy sin a 
Cy = Cjy sin a-Cy eos a 


( 2 . 10 ) 


The relationship between the two sets of aerodynamie eoeffieients is shown in 
Figure 19, where Cn = -Cz and Ca = -Cx and the side-slip is zero. The Missile DATCOM 
software determines the lift and drag eoeffieients automatieally and outputs Cl and Cd 
values direetly for the desired angles of attaek. 
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Figure 19. Relationship Between Aerodynamic Coefficients with = 0 

Source [18]: Bollino, K. P., 2006, “High-Fidelity Real-Time Trajectory Optimization for 

Reusable Launch Vehicles.” 

2. Range of Angle of Attack (AoA) 

While the fidelity of the coefficient data will be discussed further in Chapter V, an 
issue that was discovered with regard to selecting the range of AoA selected will be 
discussed here. The resulting products in subsections a and b assume a degree curve 
fit of the aerodynamic data to provide an illustration for the discussion. 

a. Initial AoA Range 

One of the inputs the user can enter into the Missile DATCOM software is the 
number and range of angle of attack values that they want the software analysis to be run 
over. When the Missile DATCOM software was first run, it was done over a range 
of -2 to 28 degrees angle of attack. The data that is outputted produces the coefficients for 
only the specific angles of attack. For optimization, an equation where the coefficients 
can be determined for any given angle of attack is needed. The desired equation was 
produced by plotting the data points and applying different curve fits to the data. 
Figures 20 and 21 show the coefficients of lift and drag and their curve fits assuming 
a 2"^^ degree polynomial fit. 
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The value shown in the plots represents the correlation coefficient. This value 
describes how well the equation represents the data. If the equation fit the data perfectly, 
the R value would be identically 1. 

The trendlines that were applied to the data produced the following normalized 
coefficient equations. 


C, =2.0809aV 1.1477a-0.0106 

( 2 . 11 ) 

= 4.2284a"-0.2776a + 0.1281 ^ 2 ) 

These equations were originally used to obtain an optimal control solution over an 
extended range of a (i.e., beyond the -2<a<28 range), shown in Figures 20 and 21. 
Note these plots only show data over a limited of attack range. 

The curve fits for -2<a<28, shown in Figures 20 and 21. This range was 
selected because it encompasses the expected autopilot limits. These curves are not 
symmetric as in Figure 18, because of the limited range in angle of attack selected for the 
data fits. 
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Figure 20. Coefficient of Lift with a 2nd Degree Polynomial Fit 

for -2 < a < 28 versus AoA 
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for -2 < a < 28 versus AoA 
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b. Extended AoA Range 

When Equation (2.11) was plotted against the entire range of angle of attacks used 
as bounds in the optimal control problem, the resulting plot (see Figure 22) showed 
increasing lift for negative angles of attack, which defies the laws of physics. 
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Figure 22. Coefficient of Lift versus Extended AoA for a 2''‘* Degree Fit 

To correct this issue, the Missile DATCOM program was run again with an 
increased angle of attack range going from -28 to 28 degrees, to include a range of 
negative coefficients within the bounds used in the optimal control problem. Figures 23 
and 24 show the new lift and drag plots with a new 2"‘^ degree polynomial curve fit over 
the increased angle of attack range. 
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Normalized Coefficient of Drag og’ Normalized Coefficient of Lift 
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tre 23. Coefficient of Lift with a Degree Polynomial Fit for 

-28 < « < 28 versus AoA 
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Figure 24. Coefficient of Drag with a 2"^^ Degree Polynomial Fit for 

-28 < cr < 28 versus AoA 
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The trendlines that were applied to the data produeed the following normalized 
eoeffieient equations. 


Q =0.0723«"+1.8962«-0.0117 


(2.13) 


= 3.6638a" + 0.0012a + 0.1108 


(2.14) 


These updated equations were used for subsequent runs of the optimal eontrol software. 


The trendlines produeed with the angle of attaek range -2 < a < 28 , Equations 
(2.11) and (2.12), were different from the trendlines produeed for the angle of attaek 
range “28 < a < 28 ^ Equations (2.13) and (2.14). The two sets of equations were plotted 
against the full range of angles of attaek and were eompared. Eigure 25 shows a 
eomparison of the two equations for eoeffieient of lift. The blue line represents 
Equation (2.11) and the red line represents Equation (2.13). Eigure 26 shows a 
eomparison of the two equations for eoeffieient of drag. The blue line represents 
Equation (2.12) and the red line represents Equation (2.14). As ean be seen, the problem 
of inereasing lift for negative AoA has been eorreeted by extending the range of angle of 
attaek data. The approximation of the drag eoeffieient is similar to before but with 
enhaneed aeeuraey for negative AoA values. 
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Figure 26. Comparison of Cd Equations 
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3. Issues on Fitting the Aerodynamic Data 

Depending on the fidelity of the model that is being used, different curve fits may 
represent the data better than others. The curve fits are an estimate of the coefficients and 
are designed to be used for simulation since tables will later be developed from testing on 
the actual missile. Thus, a lower fidelity curve fit may be an appropriate choice for initial 
simulations and optimal control runs. 

a. Coefficient of Lift 

A linear, a degree polynomial, and a degree polynomial curve fit were 
applied to the extended coefficient of lift data. They are shown in Figures 27-29. 
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Figure 27. Normalized Coefficient of Lift with a Linear Curve Fit versus AoA 
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Figure 29. Normalized Coefficient of Lift with a Degree Polynomial 

Curve Fit versus AoA 
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b. Normalized Equations for the Coefficient of Lift 
The resulting equations for the eoeffieient of lift are: 


= 1.89565a-0.0079 

(2.15) 

C, =0.0723a"+1.8962a-0.0117 

(2.16) 

C, = 2.9758a^ +0.0407a^ + 1.4393a -0.0066 

(2.17) 


Any of equations (2.15) to (2.17) eould be used to represent the eoeffieient of lift. 
In the remainder of this thesis, Equation (2.16) was used to represent the lift eoeffieient. 

c. Coefficient of Drag 

The same set of eurve fits was applied to the eoeffieient of drag plots. Beeause the 
eoeffieient of drag is parabolie, a linear eurve fit is a very poor fit for the eoeffieient of 
drag as shown in Figure 30. The 2"‘^ and 3“^ degree polynomial eurve fits are better and 
these are shown in Figures 31 and 32. As is seen, the 2"‘* and 3"^‘* polynomials are about 
the same, so either eould be used to represent the eoeffieient of drag. 
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Figure 30. Normalized Coeffieient of Drag with a Finear Curve Fit 

versus AoA 
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Figure 31. Normalized Coefficient of Drag with a 2'^'^ Degree Polynomial 

Curve Fit versus AoA 
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Figure 32. Normalized Coefficient of Drag with a 3’^‘* Degree Polynomial 

Curve Fit versus AoA 
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d. Normalized Equations for the Coefficient of Drag 

The resulting equations for the eoeffieient of drag are shown in the following 
equations. 

=-0.025« +0.3025 (2.18) 

= 3.6638a" + 0.0012a+ 0.1108 (2.19) 

Q =-0.0334a" + 3.6642a" + 0.0063a+ 0.1107 (2.20) 

In the remainder of this thesis, Equation (2.19) was used to represent the drag 
eoeffieient. Sinee Equation (2.18) does not do a good job of approximating the parabolie 
nature of the drag data, it should not be used. 

C. SUMMARY 

This Chapter introdueed the configuration and performance characteristics of the 
guided rocket model that is used in this thesis. The 6DOF and 3DOF models that can be 
used to represent the missile were introduced, as well as the assumptions for the models. 
The Missile DATCOM software used to derive the coefficients of lift and drag was 
described and the development of the coefficient of lift and drag equations was shown. 
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III. CONVENTIONAL APPROACH FOR GUIDANCE 


In order to determine the benefits of applying optimal eontrol to the missile 
guidanee problem, a baseline that the results of the optimal eontrol solution ean be 
eompared to is needed. This baseline is determined using the eonventional approaeh for 
missile guidanee. This ehapter introduees the standard guidanee law, namely proportional 
navigation [7], that is eurrently used in guided missiles. 

A. FUNDAMENTALS OF PROPORTIONAL NAVIGATION 

The fundamentals of taetieal missile guidanee are deseribed in detail in [7]. This 
seetion provides a summary of the ideas presented in [7] that are needed in order to have 
a basis of understanding of the eonventional approaeh. Proportional navigation utilizes 
the line-of-sight between the missile and the target. The guidanee law “issues 
aeeeleration eommands, perpendieular to the instantaneous missile-target line-of-sight, 
whieh are proportional to the line-of-sight rate and elosing veloeity” [7]. Mathematieally 
this ean be deseribed as: 

n^^N'V^X, (3.1) 

where ric is the aeeeleration eommand, N' is the effeetive navigation ratio, U is the 
missile to target elosing veloeity and X is the line-of-sight angle. A diagram of the 
geometry in a two-dimensional plane is shown in Figure 33. In the Figure, Rjm is the 
range from the target to the missile, Vm is the veloeity of the missile, and rit is the target 
aeeeleration. 
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Figure 33. Two-dimensional Missile-Target Engagement Geometry 

From [7]: Zarchan, P., 2007, “Tactical and Strategic Missile Guidance,” American 
Institute of Aeronautics and Astronautics, Inc., Reston, VA. 


From the geometry shown in Figure 33, the following relationships ean be 
derived [7]. The veloeity of the target ean be split into two eomponents based on the 
flight path angle of the target, fi. The components become 


V-j.^ = -Vj cos P 

'^Tz - P 


(3.2) 


The differential equations for the missile velocity and position from [7] are 

“ \yMx •> ^Mz ] 

~ [^Mx ’^Mz] 


(3.3) 


For the problem that is being examined in this thesis, the target is stationary and 
the maneuver is being conducted in the x-z plane (for maximum range). The separation 
between the missile and the target can therefore be found by 




(3.4) 
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where Xt is the position of the target and x is the eurrent position of the missile. The same 
applies in the z-axis. Given the range from the target in both of the axes, the line-of-sight 
angle ean then be found by simple geometry 

^ = tani^. (3.5) 


If the target is moving, the relative veloeity between the missile and the target is 
found by taking the differenee between the veloeity eomponents in eaeh axis. 

Vt.-Vm,] (3.6) 

Sinee the target is stationary in this problem. Equation (3.6) reduees to 

rm=[-rM, (3.7) 


From the relative veloeities and ranges, the line-of-sight rate ean be determined by 

R V -R V 

VjWx TMz TMz 


x=- 


R 


(3.8) 


TM 


where the range from the target to the missile is found by applying the distance formula 

RtM = ^IRtMx + RtMz^ • (3-9) 


Now that the range to the missile is known, the elosing veloeity ean be 
determined as [7] 


V =-R = 


{^TMx^TMx ^TMz^TMz ) 


R-r 


(3.10) 


Equation (3.10) provides the last variable that is needed to ealeulate the missile guidanee 
eommand. Onee the magnitude of the missile guidanee eommand is known, the 
aeeeleration in eaeh axis ean be determined by 


aMx=-f^c sin^ 

^Mz =^c 


(3.11) 
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B. MAXIMUM RANGE USING PROPORTIONAU NAVIGATION 


In order to implement the equations listed in the seetion above, a SIMULINK 
model was ereated for the eonventional model, shown in Figure 34. The first MATLAB 
funetion shown in Figure 34, titled PN, eontains the missile guidanee equations deseribed 
in Seetion A of this Chapter. The resulting aeeeleration that is used to adjust the missile 
LOS is outputted. The seeond MATLAB funetion embedded in the SIMULINK model 
eontains the dynamie equations for the 3DOF aerodynamie model. The aeeeleration that 
was eommanded is eompared to the aeeeleration is outputted from the missile_3DOF 
funetion bloek. The eurrent states are inputted and an updated state veetor is outputted 
along with the aeeeleration of the missile. The new position is used to determine the new 
line-of-sight angle and the proeess (eontrol loop) eontinues until the missile impaets the 
ground stopping the simulation. 

Due to the faet the problem assumes that the missile is already in flight and that 
the maneuver begins at the point when thrust turns off, the initial eonditions had to be 
developed and set to mateh the flight profile. It was assumed that the missile would be 
flying in trimmed flight at the moment the maneuver begins. Appropriate initial 
eonditions are given in Table 8. 


Table 8. Initial Conditions for Conventional Maneuver 



Units 

Value 

Xo 

m 

0 

Yo 

m 

0 

Zo 

m 

661.1 

Vo 

m/s 

514.5 

Oo 

rad 

0 

Yo 

rad 

0 


rad 

0.349 

^0 

rad 

0 
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Figure 34. SIMULINK Block Diagram for Simulating Proportional Guidance 
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In order to find the maximum range of the missile using the SIMULINK model, 
the location of the target was adjusted further and further out until the missile could no 
longer hit the target and impacted the ground shy of the intended point. The maximum 
range of the missile using the conventional approach was determined to be 8637 m with a 
flight time of 32.3 seconds. The trajectory of the missile is shown in Figure 35. Due to 
the fact that the target is stationary, the trajectory that the missile follows is essentially a 
straight line. The velocity profile versus time of flight is shown in Figure 36. 



Figure 35. Trajectory of Missile for Conventional Guidance 
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Figure 36. Missile Velocity for Conventional Guidance 


Since the boost phase of the missile is completed prior to the start of the 
maneuver, gravity is the main force acting on the missile and it assumes a ballistic 
trajectory to intercept the target. The target trajectory for a ballistic missile is 
approximately a straight line as proven in reference [7]. Due to the nature of the ballistic 
trajectory, the line-of-sight angle for the missile remain s constant for the majority of the 
flight. When missile is in close proximity to the target, the line-of-sight angle begins to 
change as the proportional navigation law adjusts the acceleration to zero line-of-sight 
error. This is shown in Figure 37. These results will be used later in the thesis to provide 
a baseline for comparison when the missile is flown using an optimal control solution. 
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Figure 37. Line of Sight versus Time 


C. SUMMARY 

This chapter introduced the conventional proportional navigation approach that 
has long been used for guidance of rockets. In addition, a maximum range maneuver was 
simulated to establish a baseline performance that will be used to compare against the 
results achieved through optimal control. 


50 




IV. OPTIMAL CONTROL THEORY 


There are multiple paths that a missile ean take to reaeh the desired end state. 
Eaeh path will require a different amount of resourees, sueh as fuel or time. Optimal 
eontrol theory looks to find the path that best minimizes or maximizes a given objeetive 
funetion with respeet to the eonstraint set on the problem. For example, a user may want 
the quiekest route. The objeetive funetion in this ease is to minimize time of flight. 
Conversely, a user may want to extend the range of a missile so the objeetive is to 
maximize the distanee travelled form the launeh point. In either ease, the optimal eontrol 
solution will yield the path that optimizes the objeetive. 

This ehapter provides the basie baekground needed for this thesis in regard to 
optimal eontrol theory. A more thorough diseussion of optimal eontrol ean be found in 
[15]. Referenees [18], [19], [21] also provide a good introduetion to the basies of optimal 
eontrol. This ehapter eontains a summary of the important eoneepts from these 
referenees. In order to solve optimal eontrol problems and find the appropriate path for a 
given task, an optimal eontrol software tool is typieally used. This is beeause solving the 
problem analytieally is often intraetable. For this thesis, DIDO was used to obtain all the 
optimal eontrol solutions. This ehapter, therefore, also ineludes an introduetion to the 
DIDO psuedospeetral optimal eontrol software. 

A. THE OPTIMAL CONTROL PROBLEM 

The generie optimal eontrol problem as derived in [15], [18], is shown in 
Equation (5.1). Within this equation, J is the eost funetional, E is the endpoint eost, F is 
the running eost. The value of the eost funetional is subjeet to a set of state variables, x, 
and eontrol variables, u. The equation x = f(x,u), in Equation (5.1), defines a set of 
dynamie equations for the problem. In addition e defines the set of end point eonstraints 
(boundary eonditions) for the problem. 
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iPo) 


! 

MinimizeJ =E{Xj,tf) + ^F{x{t),u{t))dt 


Subject to ; 

x{t) = f{x{t),u{t)) 

x{tj = 


(4.1) 


t=f 


t, =E 


e{xit.)) = 0 


Pontryagin’s Principle is stated as follows: “Given an optimal solution to Problem 
Po, there exists a eoveetor fimetion, A(»), and a eoveetor, v, that satisfy the Hamiltonian 
minimization eondition, the Hamiltonian value eondition, the Hamiltonian evolution 
equations, the adjoint equations, and the transversality eondition” [15]. These eonditions 
lead to the formulation of a boundary value problem (BVP), whieh ean be solved to find 
the optimal solution. 

The first step in this proeess is to develop the Hamiltonian funetion which is 
defined as: 

H{A,x,u) = F{x,u) +f{x,u) t4 2i 


The Hamiltonian introduces the eostates A , whieh represent the adjoint eoveetor. The 
eoveetor has the same dimension as the state veetor, x, but does not oeeupy the same 
spaee asx [15]. The adjoint dynamies are defined as: 

= (4.3) 

OX 

Aceording to Pontryagin’s principle for the control function u(») to be optimal, 

u(») must globally minimize the Hamiltonian for every t from to to (/, while holding 

A and x eonstant at each instant in time [15]. This is known as the Hamiltonian 
Minimization Condition (HMC), whieh is defined as: 
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Minimize 

{HMC)\ 

[Subject to 


ueU 


(4.4) 


If the controls do not have any constraints, the HMC is reduced to the Euler- 
Lagrange equation which is given as; 



du 


(4.5) 


In the presence of constraints on u, partial derivative 


dH 

du 


must be interpreted as a 


switching function, that is: 


5 


u 


dH 

du 


(4.6) 


The sign of then determines the value of the control variable at any instant in time. 


The optimal control BVP should have 2N boundary conditions where N is the 
number of state variables. In order to find any missing boundary conditions, the Terminal 
Transversality Condition given in Equation (4.6) is used [15]. 



dE 

dx^ 


(4.7) 


where E is the endpoint Eagrangian, defined as 


E{v, = E(x(tj)) + o^e(x(tf)) 


(4.8) 


While the Transversality condition provides Nx new equations, it also produces Ne 
unknowns. There is one final missing boundary condition which comes from the 
Hamiltonian Value Condition, given in Equation (5.8). 

Equations (4.7) and (4.9) provide the remaining needed boundary conditions in 


order to be able to construct and solve the boundary value problem (BVP) associated with 
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the optimal control problem. A challenge in solving the BVP is that it cannot typically be 
done analytically. Thus, a software tool is needed. 


B. INTRODUCTION TO DIDO 

DIDO is a MATLAB-based software package that was originally developed by 
Fahroo and Ross of the Naval Postgraduate School to solve complex optimal control 

problems [22]. DIDO was “originally developed to make it “easier to solve optimal 
control problems” by allowing the user to focus on defining the problem they want to 
solve as opposed to worrying about how to solve it” [23]. This tool takes the system 
dynamic equations and boundary conditions of the basic optimal control problem as 
described in Section A and develops a solution through use of pseudospectral optimal 
control theory. A more thorough discussion of the theories and mathematics behind 
DIDO can be found in [16]. 

DIDO requires the user to state the problem in a specific format that is encoded in 
four MATLAB m-files. The four function files are the “dynamics,” “cost,” “events,” and 
“path”-files. The cost, dynamics, and event files are mandatory for every problem but the 
path file is optional. The cost function describes the endpoint cost and running cost for 
the problem. The cost for the problem is determined by the user and is set based on what 
the problem is intended to minimize (e.g., time) or maximize (e.g., range). The dynamics 
function contains the differential equations or state-space model to be optimized. In the 
event function m-file, all the known boundary conditions for the problem are entered. 
The path function allows the user to provide path constraints at the nodes. The path is an 
optional file and is only used if the problem has path constraints. The problem being 
examined in this thesis does not require use of the path function. However, a path 
constraint could be used if it is necessary to limit some aspect of the missile performance 
so that it is consistent with operational requirements, for example, bounding the 
maximum acceleration. 

The m-file that is needed in order to be able to call DIDO in MATLAB is a main 
problem file. Within the main file, the upper and lower limits for the state variables, 
control variables, and clock times are defined. In addition, the number of nodes that will 
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be used to determine the solution (the eomputational grid) are also defined in the main 
file. The number of nodes determines the aeeuraey of the solution. One advantage to 
DIDO is that it ean be used iteratively and a lower fidelity solution ean be used and as a 
seed to refine the solution [21]. This provides a way to improve the aeeuraey of the 
solution quiekly. The main file also allows the user to output any desired data or plots 
from the run. Onee the eost, dynamies, events, and main files have been eompleted, 
DIDO ean be run. DIDO performs analysis of the problem and if a solution is found, 
DIDO will output a message indieating that an extremal solution has been found. A 
solution is extremal if it is a minimizer for the HMC [22]. While Pontryagin’s boundary 
value problem is not direetly eoded, DIDO provides data on the eostates that ean be used 
to eheek of the solution against Pontryagin’s neeessary eonditions as part of a verifieation 
and validation (V&V) step. 

For any extremal solution, it is important to verily and validate the results that 
were produeed. This ensures that an optimal solution was in faet obtained and that the 
problem that was solved was indeed the one intended by the user. In addition to eheeking 
Pontryagin’s neeessary eondition, the feasibility of the results should be eheeked by 
utilizing the built in MATLAB ODE tools. In this V&V test, the outputs from the DIDO 
run are propagated through the ODE funetion (dynamies) and the results eompared in 
terms of their residual differenees. An easy way to do this is by ereating a simple 
eomparison graph where both the outputs from the DIDO run and the propagation are 
plotted together. When the results are eompared, if the optimal eontrol results matehes 
the results of the ODE propagation within some aeeeptable toleranee, the solution ean be 
eonsidered feasible. 

1. Scaling the Problem 

In any numerieal algorithm, the sealing of the problem is also important. If the 
problem is badly sealed, the algorithm will be unable to solve the problem. This is 
typieally seen in DIDO when the relative magnitudes of the states and the eostates are 
large. This problem ean be addressed by refining the solution through sealing. The 
sealing is aeeomplished through the use of “designer units” [23]. Equation (2.7) is the 
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dynamics equations for the 3DOF model. The dynamie model is formulated in state spaee 
and is given byi: = f{x,u,t). In this state spaee, the eontrol funetion u is the input into 
the dynamies equation to produee a response. The states for the model ean be assembled 
into a veetor as; 


=[x,y,z,v,(j,y] 


(4.10) 


The basie idea behind designer units is that eaeh element in the state veetor ean be sealed 
by an arbitrary number seleeted by the user as shown in Equation (4.11), 


—T 

X 


X y z V (J y 

x’t’z’f’e’t 


whereX, Y, Z, ete., are the designer seale faetors. 


(4.11) 


The same method of sealing applies to the eontrol variable and the time units, 
whieh beeome; 


a jd 

a’m 


(4.12) 



(4.13) 


While the sealing faetors are arbitrary, the goal is to piek “good” designer units. 
“A fundamental rule for sealing optimal eontrol problems is to ehoose designer units in 
sueh a manner that the sealed states and eostates are roughly the same order of 
magnitude” [23]. The reason it is important for the states and eostates to be roughly the 
same order of magnitude is due to the relationship between the two. 


state UNITS 


(4.14) 


For the 3DOF missile problem the sealed eostate veetor beeomes; 

X rj-i ‘ y rp ‘ Z p ‘ V p ‘ (7 p ‘ 7 p 


(4.15) 
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Pontryagin’s Principle and the Hamiltonian Value Condition do not ehange when 
the Hamiltonian is sealed whieh is why problem sealing is effeetive. For a praetieal 
problem, the sealing is applied within the MATLAB eoding and a trial and error 
approaeh is used in order to find the values that produce a solution where and the states 
and eostates at approximately the same order of magnitude [23]. 


2. Coding the Scaling in MATLAB 

In order to apply sealing to the problem, several additional lines of code are 
required in the m-files that are to be used by DIDO. The first step is to create a global 
variable by the name of SCALES in the main file. This is done for convenience so that 
the scale factors can be easily shared between the different m-liles. In the main file, the 
scales are then defined for each of the states, controls, and time variable as follows: 


global SCALES; 
SCALES.X= 1; 
SCALES.Y= 1; 
SCALES.Z= 1; 
SCALES.V= 1; 
SCALES.SIGMA= 1; 
SCALES.GAMMA = 1; 
SCALES.ALPHA= 1; 
SCALES.MU= 1; 
SCALES.T= 1; 


(4.16) 


In the dynamies file, in the preamble where the states, eontrols and nodes are 
defined from the primal data strueture, the sealing faetor for eaeh variable is applied. 
Then when defining xdot, the veetor is multiplied by the time seale faetor and then eaeh 
element of the veetor is divided by its sealing faetor returning the output to the original 
value. Thus, the dynamies are computed in engineering units but are operated on by the 
numerieal algorithm in terms of designer units. In the main file, all of the bounds for the 
states, controls, time and events must also be divided by the seale for the appropriate 
value. Of note, the scaling used when running the program with lower nodes may need to 
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be refined as the results are stepped into higher node runs for increased accuracy. An 
additional benefit to the scaling is that it also allows for quicker run times. 

As an example, for the optimal control problem that is solved in later Chapter V, 
when DIDO is run without any scaling the output for the states x and y versus time are 
shown in Figures 38 and 39. As shown, the ranges of the two states differ by five orders 
of magnitude so scaling can be used to collapse this range. 



Figure 38. Down-Range Distance versus Time 
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Figure 39. Flight Path Elevation Angle versus Time 

The scaling not only applies to the relationship between the states, but also the 
relationship between the state and it’s costate. The costates for x and y versus time are 
shown in Figures 40 and 41. The range of x and differ by four orders of magnitude 
and the range of y and differ by three orders of magnitude. The ideal range between 

the states and costates is a 1:1 relationship with the range of both values being of the 
same order of magnitude. 
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Figure 40. versus Time for Maximum Range Maneuver 



Figure 41. versus Time for Maximum Range Maneuver 


Designer units can be applied to these variables in the dynamics files as 


X = primal.states(l,:)*SCALES.X; 
sigma = primal.states(5,:)*SCALES.SIGMA; 


(4.17) 
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This same principal is applied to all of the other states. Since the dynamics m-file must 
output the scaled dynamics {xdof) for the algorithm to operate on, the dynamics 
(calculated in engineering units) must be converted back to designer units. Converting the 
dynamics from engineering units to designer units is done as shown in Equation (4.18). 


xdot(:,i) = SCALES.T*[xDot(i)./SCALES.X; 

sigmaDot(i)/SCAEES.SIGMA]; 


(4.18) 


When DIDO is run with the scaling applied, the resulting scaled states and costates are 
shown in Figure 42. Referring to Figure 42, the scaled states and costates are now similar 
in order of magnitude, which improves the quality of the numerical solution. 



Figure 42. Scaled States versus Time for Maximum Range Maneuver 


C. SUMMARY 

This chapter provided the basic background needed for this thesis in regard to 
optimal control theory. This chapter also introduced the DIDO psuedospectral optimal 
control software, the software package used to solve the optimal missile trajectories in 
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this work. Issues related to numerieal eonditioning via sealing through designer units 
were also diseussed. The optimal eontrol framework will be applied to the missile 
problem in the next ehapter. 
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V. MAXIMIZING THE RANGE OE A GUIDED ROCKET 


In this chapter optimal control theory is used to obtain a guidance solution to 
maximize the range of the fictitious missile under study. The bounds and conditions that 
were used for the problem setup as well as development of the necessary conditions for 
optimality will be shown. The results of the solutions obtained using DIDO will be 
presented and discussed. The resulting maneuver will be compared to the conventional 
proportional navigation based solution. As will be seen, the application of optimal control 
theory allows the range of the missile to be enhanced significantly. 

A. DESCRIPTION OF THE PROBLEM 

The desired outcome of the optimal control problem studied in this chapter is to 
maximize the range that the missile is able to travel. In order to start building the optimal 
control problem formulation, the 3DOF equations developed in Chapter II will be used. 
As mentioned previously, for the scope of this thesis, the problem will consider a guided 
rocket traveling in the mid-to-terminal phase of flight only. Thus, it will be assumed that 
the problem begins once the rocket motor is burned. The 3DOF optimal control problem 
was set up utilizing the version of the dynamics equations shown in Equation (2.7). The 
state and control variables for the problem are repeated in Table 9 for convenience. 
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Table 9. Reduced 3D0F Model States and Controls 


State 

Symbol 

Units 

x-position (down-range distance) 

X 

m 

y-position (cross-range distance) 

y 

m 

z-position (altitude) 

z 

m 

velocity (airspeed) 

V 

m/s 

flight path heading angle 

a 

rad 

flight path elevation angle 

Y 

rad 

Control 

Symbol 

Units 

angle of attack 

a 

rad 

bank angle 

P 

rad 


Since the missile’s autopilot utilizes inner loop control logic, the control variables 
for the optimal control problem can be written as 

d=u„ 

(5.1) 

which allows the optimal control solution to accommodate the time constraints from the 
inner loop. 

The maximum range problem is solved by minimizing the negative value of the 
final downrange position, . Hence for the maximum range problem, there is an 

endpoint cost but no running cost. Problems with a specified crossrange are not 
considered. The maximum range optimal control problem is therefore stated as 
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V = -g sin / 

m 

L sin // 
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mv cos / 
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d=u^ 
li = u^ 


= -X, 


t^=t 


0 


e(x(t^)) = [z^_ 


u <\u \<u 

\ a \ a 


(5.2) 


In Equation (5.2), jc° is a given vector that consists of the initial conditions for the 
states at time to. The values that comprise are given in the next section in Table 10. 

Equation (5.2) also states that e(x(t^)) = [z^ J (i.e., the maximum range is reached when 

the missile impacts the ground). The numerical values that comprise the initial and final 
conditions are given in the next section in Tables 11 and 12, respectively. Since the 
problem statement does not include a running cost, the value is equal to zero at the final 
time. 


B. APPLICATION OF THE OPTIMAL CONTROL THEORY 

Now that the problem has been stated, the next step is to apply the process that 
was described in Section A of Chapter IV, in order to derive the necessary conditions for 
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optimality. The first step is the formulation of the Hamiltonian utilizing Equation (4.2). 
The Hamiltonian for the Problem (5.2) is 


H[X,x,u) = [veosyeoscr] + /l^ [veos/sincr] + /l^ [vsiny] 


+x 


\-D . ] 

+ 4 

L sin// 

-gsm/ 

m J 


mv cos/ 


(5.3) 


L g 

—eos//-—cos/ 
mv V 


+ 


Taking the partial derivative of the Hamiltonian with respect to the control 
variables gives 


dff 

du^ 


= X=Q 


dH 

du^ 



(5.4) 


Thus, the Hamiltonian Minimization Theorem must be applied, where 


I Minimize 
{HMC)\ 

[Subject to 


H{X,x,u) 

<u<u^ 


(5.5) 


due to the fact that the control variables appear linearly in Equation (5.3). The equations 
in Equation (5.4) become switching functions, as shown in [19], where 


5 


a 


4>o 

, 4=0 






<U„ < U„ 


(5.6) 


The same principal can be applied to the other control variable, which produces 
the second switching function 
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(5.7) 




4<o 

^=0 


>“min 


U —U 

fi /^max 


^ < M. 

-“max 


These switching functions given in Equations (5.6) and (5.7) can be utilized as an 
additional means for validation and verification of the optimal control solution. For 

example, and can be plotted to check that if is at the upper bound, then ^ ^ , 

etc. The same analysis can be performed with and u^. If the trajectories obey the 

inequality for the switching function, then the necessary condition for optimality is 
satisfied. If not, then the solution may not be correct and the reason for the discrepancy 
must be isolated and corrected. 

Because the adjoint equations are only being used for verification and validation 
rather than problem solving, p can be assumed constant in the drag and lift equations 
when the costate dynamics are derived. The application of Equation (4.2), the adjoint 
equation becomes 
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(5.7) 


From Equation (5.7), since A^=A^=A^=0, the eostates for the x, y, and z 

variables must be eonstant. This provides another pieee of information that ean be used to 
verify the DIDO outputs. If the numerieal values of the eostates are not eonstant, then the 
solution ean be determined to be invalid. 

The transversality eondition is applied next. The Endpoint Eagrangian beeomes 

E(v, x(tj)) = E {^Xj ) + v^e{Xf) = -x^ + , (5.8) 
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When partial derivative of the Endpoint Langrangian is taken with respeet to eaeh of the 
states (at ), the value of the eostate at time is obtained. These values, together with 

the behavior of the adjoint equations ean also be used for verifieation and validation of 
the optimal eontrol solution. 
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i(>,) 



dXf 



dy/ 


dz, 


■ = t>, 


K{‘,) 

K(‘,) 


dE 

dVf- 

dE 


= 0 


= 0 


d(jf 



dXf 


(5.9) 


1. State and Control Variable Bounds 

As mentioned in Chapter IV, Seetion B, the main program file for the DIDO 
problem eontains the required bounds for eaeh state variable. The eonstraints serve to 
provide a restrietive range to seareh for a feasible solution to the optimal eontrol problem. 
While the goal is to limit the seareh range, the seareh range needs to be large enough that 
it eneompasses the range in whieh a feasible solution ean be found. 

For the problem addressed in this thesis, the x-limit is based on having a range 
that allows the entire maneuver to be eompleted. The x-limit was adjusted until the final 
solution was inside the upper bound ensuring that the maximum range of the missile 
eould be reaehed. The z-limit set the eeiling for the maneuver. The allowed eross-range 
was ±100 m. The veloeity limits ranged from the maximum speed of the missile to 100 
m/s whieh was estimated to be the minimum speed needed on missile impaet at the target. 
The bounds for the states and eontrols are shown in Table 10. 
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Table 10. Bounds for the States and Controls for Maximum Range Maneuver 


State 

Lower Bound 

Upper Bound 

Units 

X 

0 

30000 

m 

y 

-100 

100 

m 

z 

0 

3000 

m 

V 

100 

625 

m/s 

a 

-71/6 

71/6 

rad 

Y 

-7C/6 

7C/6 

rad 

Control 




a 

-0.349 

0.489 

rad 

h 

-0.0872 

0.0872 

rad 


-0.349 

0.489 

rad/see 


-0.0872 

0.0872 

rad/see 


2. Initial Conditions 

As stated above, the maximum range maneuver begins with the missile in-flight at 
the moment that the engine has burnt out and the thrust goes to zero. The missile is 
assumed to be in steady state eontrolled flight at the moment that the terminal maneuver 
begins. The kn own initial eonditions are shown in Table 11. The initial eonditions are the 
same eonditions that were utilized for the eonventional approaeh in order to be able to 
properly eompare to the results from the OCP. 
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Table 11. Initial Conditions for Maximum Range OCP 



Units 

Value 

Xo 

m 

0 

Yo 

m 

0 

Zo 

m 

661.1 

Vo 

m/s 

514.5 

Oo 

rad 

0 

Yo 

rad 

0 

Oo 

rad 

0.349 

^0 

rad 

0 


In order to determine the initial value for the angle of attaek, the faet that the rate 
of change of the flight path elevation angle would need to be zero for the missile to be in 
straight and level flight was used. From Equation (2.7) it is given that 

;k = -^cos//-— cos/ (5.10) 

mv V 

Substituting the known values at the initial condition, the equation reduces to 

--g = 0, (5.11) 

m 

where 

(5.12) 

By substituting in the equation for Cl that was developed earlier in Chapter II, the 
angle of attack for straight and level flight was determined to be 0.349 rad (20 deg). For 
the optimal control model, the missile is being controlled through adjustment of the angle 
of attack. On the physical missile, however these changes are implemented through 
applying fin deflections to produce the desired results. Thus, for implementation of the 
optimal maneuver, an inner loop autopilot would need to be used. 
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3. 


Final Conditions 


In addition to the initial conditions, the events function allows the user to input 
any known final or endpoint conditions. For this problem, the distance that the missile 
would travel was unknown as well as the velocity at the point of impact. In addition, the 
orientation of the missile with respect to the reference plain at the time of impact was not 
known. The only known endpoint condition was that missile would impact the ground, 
thus z/= 0 as listed in Table 12. 

Table 12. Known End Point Conditions for the Maximum Range OCP 



Units 

Value 

Zf 

m 

0 


C. MAXIMUM RANGE SOLUTION 

With all the data defined for the cost, event, dynamics, and main files, the 
problem case was coded and solved using DIDO. 

I. Selection of Coefficient of Lift and Drag Fidelity 

In Chapter II, equations for the coefficient of lift and drag were determined from 
the Missile DATCOM. Equations (2.15)-(2.20) provided three different curve fits that 
could be used for each of the coefficients. When DIDO was initially run, a third degree 
polynomial curve fit was used because it had the highest correlation value of the three 
fits. The results of these initial runs indicated that the problem was infeasible. 

While the results from a DIDO run for an infeasible solution cannot be trusted as 
accurate, one observation that was noticed was that the costates were incredibly large, 
shown in Eigure 43. This indicated that problem scaling was likely the reason why the 
problem was being determined infeasible. 
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Figure 43. Costates versus Time for DIDO Run with a Degree Polynomial 

Curve Fit Applied to Cl and Cd 

The two eostates that were drastieally larger than the other costates were the one 
for velocity and flight path angle. When examining the equations for the costates, 
Equation (5.5), it was noted that they were the only two costates that contained the 
coefficients of lift and drag. In addition, both equations had cosy orcos^y in the 
denominator. For small angles of attack, the result of the angle of attack being cubed in 
the equation caused the value for gamma to approach nil, which produced a very small 
number in the denominator causing the values of the costates to become quite large. 

To test the theory, the problem was run with a linear curve fit for the coefficient 
of lift and drag equations and DIDO produced a feasible result with the costates all within 
the same order of magnitude. The problem was then run again with a second degree curve 
fit and DIDO once again produced a feasible result with the costates within the same 
order of magnitude. This process provided a valuable lesson learned. When utilizing a 
model, such as Missile DATCOM output, to produce equations for estimating the 
coefficients of lift and drag it is best to start with a lower fidelity model for initial 
simulation runs. Once the problem has been properly scaled, the fidelity can then be 
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increased incrementally the sealing adjusted appropriately at each increment. For the 
results that will be discussed later in the chapter, both the coeffieient of lift and drag 
equations are the equations obtained for seeond degree polynomial curve fits. 

2. Results 

a. Scaling the Problem 

When the problem was run with the order curve fits for the aerodynamic data 
(Equations 2.16, 2.19), a feasible solution could be obtained using A^ = 16 eolloeation 
nodes. However, the resulting eostate for the flight path angle. Figure 44, was still several 
orders of magnitude greater than the other eostates. In addition, the variation in the 
Hamiltonian, Figure 45, was large, whieh is ineonsistent with the assoeiated necessary 
eondition. As stated in Chapter IV, Seetion B, the way to resolve these issues is through 
proper scaling. 
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Figure 44. Costates versus Time for Unsealed Maximum Range Maneuver 
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Figure 45. Flamiltonian for the Unsealed Maximum Range Maneuver 

In order to determine the proper “designer units” for this problem, multiple trial 
run were eondueted. The resulting sealing faetors for the problem were determined as 
follows: 


SCALES.X= 1000; 

SCALES.Y = 50; 

SCAEES.Z= 100; 

SCAEES.V = 50; 

SCAEES.SIGMA= 1; (5.13) 

SCAEES.GAMMA = 0.1; 

SCAEES. ALPHA = 1; 

SCALES.MU= 1; 

SCALES.T=1; 

With these designer units the ranges of the states and costates became comparable (as 
will be shown next) and the problem run time improved. 
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b. Solutions for the Scaled Maximum Range Problem 

The problem was then re-run and refined using = 16, 24, and 32 collocation 
nodes utilizing the scaling factors developed in the previous section. After DIDO 
produced a feasible solution, the states, costates, controls, Hamiltonian, and cost for the 
problem were outputted. The position of the missile versus time is shown in Figure 46. 
Since the cross-range is zero (even though this condition was not specified), the results 
show that the range is maximized by flying in the x-z plane. 



Figure 46. Position versus Time for the Maximum Range Solution 

When the maneuver was initiated, the missile is seen to initially climb to an 
altitude of 1954 m before beginning its decent. The missile impacted the ground at a time 
of 94.02 seconds. The missile’s altitude throughout the flight is shown more clearly in 
Figure 47. As mentioned above, the optimal trajectory for the missile was in the x-z- 
plane. The trajectory of the missile in this plane is shown in Figure 48. 
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Figure 47. Altitude Profile for the Maximum Range Maneuver 
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Figure 48. Trajeetory of the Missile for Maximum Range Maneuver 
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The missile was traveling at a speed of 121.4 m/s when it impacted the ground, 
as shown in Figure 49, which exceeds the minimum required velocity of greater than 
100 m/s for impact. 
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Figure 49. Velocity versus Time for Maximum Range Maneuver 


The angle of attack for the missile was set with an upper limit of 0.489 rad 
(28 deg). The missile reaches the upper limit that was established near the end of the 
maneuver, shown in Figure 50. This behavior is similar to the conventional guidance 
approach which commands a large normal acceleration near the end of the flight in order 
to zero the line-of-sight error as shown previously in Figure 37. While the missile is 
climbing and then descending in this maneuver, the missile traveled in a zero crossrange 
path which resulted in negligible bank angles being applied to the missile, shown in 
Figure 51. If the desired maneuver required a turn or the missile was intercepting a 
moving target, this control would be utilized. 
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Figure 50. Angle of Attack versus Time for Maximum Range Maneuver 



Figure 51. Roll Angle versus Time for Maximum Range Maneuver 

The costates for the problem with scaling applied are shown in Figure 52. When 
compared to Figure 42, the effects of proper scaling on the costates can be clearly seen. 


79 





As mentioned above, Equation (5.7) showed that should be eonstant. 

Figure 52, shows this to be true and is a validation point for the optimal control solution. 
When the costates listed above are plotted against time. Figure 52, the values are not 
exactly constant (due to numerical tolerances) but are within an acceptable level of error. 
As the fidelity of the model is increased through the bootstrapping to a higher number of 
nodes, the costate becomes smoother and settles towards being a constant value. 
However, even with infinite precision there will still be some error in the model. The 
costates for the other states are shown in Figure 53. As noted earlier, the behavior of 
these costates is complex so these trajectories are not useful as far as V&V is concerned. 
However, trasversality provides some conditions that can be checked. (Rachel pis. 
confirm that the behavior of the remaining costates (at the end time) is correct with 
respect to the transversality conditions.. For example. Equation (5.9) states 
that = , and that 2,^ j = 0 , etc. The results obtained from the optimal control 

solver are consistent with these equations. 



Figure 52. Scaled Position Costates versus Time for 
Maximum Range Maneuver 
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Figure 53. Scaled Costates for Velocity and Flight Path Angles versus 

Time for Maximum Range Maneuver 

The Hamiltonian for the scaled problem is shown in Figure 54. The variation in 
the Hamiltonian is nearly zero, as required. The fact that the solution satishes these 
necessary conditions allows one to conclude that an optimal solution has indeed been 
found and that the missile’s range has been maximized per to problem limits and 
constraints. 
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Figure 54. Hamiltonian versus Time for Maximum Range Maneuver 
c. Verification and Validation 

In order to verify that the optimal control solution is indeed a dynamically 
feasible solution, the next step was to perform a propagation test [15]. The optimal 
controls that were produced were propagated through the dynamic equations for the 
problem, Equation (2.7), using the ODE45 solver in MATEAB. The algorithm uses the 
control and time values associated with the nodes and interpolates to find the value for 
the control needed at a specific time. The results of the propagation were plotted against 
the optimal control solution for all of the states, shown in Figures 55-58. As can be seen 
the propagated solution closely matches the solution obtained from DIDO. Thus, the 
solution is feasible for implementation. 
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Figure 55. Verification of Position versus Time along the X-, Y- and 
Z-Axis for Maximum Range Maneuver 



Figure 56. Verification of Missile Flight Path for Maximum Range Maneuver 
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Figure 57. Verification of Velocity Vector for Maximum Range Maneuver 



Figure 58. Verification of Flight Angles for Maximum Range Maneuver 
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d. Impact of AoA Range 

In Chapter II, Section B, Paragraph 2, the range of AoA used to determine the 
coefficient of lift and drag plots within the MissileDATCOM was discussed. To show the 
impact that the incorrect range (and hence incorrect aerodynamic data) can have on the 
final solution, the results of the optimal control problem using Equations (2.10) and 
(2.11) are compared to the results obtained using Equations (2.12) and (2.13), shown in 
Figure 55. 



0 0.5 1 1.5 2 2.5 

Down-Range Distance (m) x 10“* 

Figure 59. Comparison of Missile Trajectory Utilizing Curve Fits Derived 

from Different AoA Ranges 


When the increased range of angles of attack to include the full range of the 
negative angles of attack, the distance the missile traveled increased significantly. 
Utilizing the incorrect solution resulted in underestimating the missile range by 3328 m. 
This result is included to be utilized as a lessoned learned for future work on this topic. 
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e. Setting Upper Bounds for AoA 

In the results shown in subseetion b, Figure 48, the angle of attaek trends toward 
the upper bounds limit toward the end of the flight. As mentioned in Section B, the upper 
bound was set to 0.489 rad (28 deg). It was assumed that to ensure penetration of the 
target, the angle of attack would need to be less than 30 degrees when impacting the 
ground. Nonetheless, the effect of the upper bounds of the angle of attack was examined 
to determine the effect these limits could have on the optimal solution. The trajectories of 
the missile utilizing four different angle of attack limits are shown in Figure 60. 



Figure 60. Trajectories of the Missile for a Maximum Range Maneuver 

with Various AoA Limits 

The trajectories were all very similar and the path the missile follows (shown in 
Figure 60) appears to be nearly the same for all four angle of attack limits. When the 
point of impact for the missile is zoomed in on the graph (see Figure 61) it is seen that 
there is a difference of 350 m between the maximum angle of attack at 0.349 rad (20 deg) 
and 0.785 rad (45 deg). Depending on the situation, the extra 350 m in range could be a 
factor that warrants the consideration of adjusting the AoA limit to a larger value. This 
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decision, however, will depend on the constraints on missile effectiveness, aspects that 
are not considered here. 



Figure 61. Trajectory of Missile in Last Portion of Flight 

While the differences in the range of the missile are rather small compared to the 
overall range in the terminal phase, the impact of the maximum angle of attack bound can 
be seen the most in the velocity at the end of the flight, shown in Figure 62. When the 
maximum angle of attack is limited to 0.349 rad, the final velocity at the point of impact 
is 159.8 m/s, while a maximum angle of attack at 0.785 rad the final velocity is at the 
lower bound at 100 m/s and the time of flight is only a few seconds shorter. 


87 










550 



Figure 62. Velocity versus Time for Varying Maximum AoA 

Depending on the target that is selected, a slightly (350 m) shorter range with a 
higher (+60 m/s) impact velocity may be desired over the increased range and slower 
impact speed. The higher velocity may allow the missile to penetrate further prior to 
detonation if the desired target is a building or bunker. If a shallower angle of attack on 
impact and decreased velocity is permissible, then the increased range may be the desired 
effect. 

The comparisons shown above were done using the coefficient of lift and drag 
equations that were developed over a range of -28 < « < 28 . In aerodynamic flight, the 
attack angle produces an increase in lift until the angle of attack has increased to a point 
where the relative surface area of the missile to produce lift has decreased enough that the 
lift for the body begins to decrease. The larger angle of attacks examined may be beyond 
this point which would result in even further decreased ranges and impact velocities. The 
effect described above was not reached within the bounds for the optimal control problem 
and determination of the angle where this begins to occur is outside the scope of this 
thesis. 
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D. COMPARISON OF CONVENTIONAL AND OPTIMAL SOLUTIONS 


1. Maximum Range Maneuver 

When the trajectories of the missile are plotted against each other for both 
methods, shown in Figure 63, the difference between the two flight paths is substantial. 
Both of the problems began with the missile in the same steady flight condition but the 
fundamentals behind each approach make a significant difference in the outcome. The 
conventional approach is driven by a simple law involving the line-of-sight and range to 
the target. In the terminal phase of flight, this results in the missile assuming a ballistic 
trajectory. The optimal control solution, on the other hand, causes the missile to initiate a 
climb in the beginning of the maneuver utilizing the initial airspeed of the missile to gain 
altitude. 



Figure 63. Comparison of the Missile Trajectory for Conventional and 
Optimal Control Method for Maximum Range Problem 

The maximum downrange distance that could be achieved by the conventional 
approach to missile guidance with the given initial conditions was 8637 m. By utilizing 
the optimal control solution, the maximum range of the missile increased to 24850 m. 
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This was an improvement in range of 16213 m (8.75nm) which is nearly a 300% increase 
in the maximum range of the missile over the conventional proportional navigation based 
approach. 

The time of flight was also significantly different from the conventional to 
optimal solution, shown in Figure 64. The time of flight for the missile using the 
conventional approach was 32.3 seconds, while the time of flight for the missile utilizing 
the optimal control solution was 94.02 seconds. This was an increase in the overall time 
of flight by a factor of 2.91. The increased range and time of flight allows the missile to 
be fired from further away and provides increased time for the aircraft firing the missile 
to get clear of the area. Both of these benefits can potentially improve the chances of 
mission success. 



Figure 64. Time of Flight Comparison for Conventional and Optimal Control 

Methods for Maximum Range Problem 

For both the conventional and optimal control approach, the velocity was plotted 
against time, shown in Figure 65. While the conventional approach has a much steeper, 
almost linear decrease in velocity, the optimal control solution has a gradual decrease 
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with the velocity leveling out slightly as the missile was at the top of the parabolic path 
and began its decent. In addition, the impact velocity was similar for both the 
conventional and optimal control solution, which means that increasing the range did not 
have an adverse impact on the terminal kinetic energy of the missile. 



Figure 65. Comparison of Velocity versus Time for 
Maximum Range Problem 


2. Minimum Time Solution 

Because the maximum range produced with the optimal control solution was 
much greater than the maximum range that could be achieved by using the conventional 
approach, the range itself was the only aspect that could be clearly compared between the 
two approaches. To provide additional comparison of the two approaches for missile 
guidance, an additional optimal control problem was solved. 

a. Problem Statement 

For this problem, the desire was to minimize the time of flight to the same range 
obtained using the conventional approach. The second optimal control problem took the 
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same initial conditions as were used in the initial problem, shown in Table 11. However, 
an additional end point condition was added as shown in Table 12, by setting the 
downrange distance to the maximum range that was achieved by the conventional 
approach. The cost of the problem was also changed to minimize the time it would take 
for the missile to impact the target. 


Table 13. Known End Point Conditions for Minimal Time OCP 



Units 

Value 

Xf 

m 

8637 

Zf 

m 

0 


The optimal control problem for the minimum maneuver is as follows: 


Minimize 


{Pc 


OC Min. Time 


)= 


= C 


Subject to: 

X = vcosycoscr 
y = vcosysincr 
z = V sin y 

. -B 

V = -g sm y 


cr = 


m 

L sin ju 
mv cos y 


L g 

y = —cos//-—cosy 
mv V 

d = u^ 

// = u^ 


to =t 


0 


e(x(ty)) = [x^ -8637,z^ J 


/ ^ I I ^ u 
u <\u \<u 

a \ a\ cc 

L I/ U 


(5.14) 
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The bounds of the states and eontrols were also adjusted for the new maneuver as 
shown in Table 13. 


Table 14. Bounds the States and Controls for Minimum Time Maneuver 


State 

Lower Bound 

Upper Bound 

Units 

X 

0 

10000 

m 

y 

-100 

100 

m 

z 

0 

2000 

m 

V 

100 

625 

m/s 

a 

-7c/6 

7c/6 

rad 

Y 

-7C/6 

71/6 

rad 

Control 




a 

-0.349 

0.489 

rad 

p 

-0.0872 

0.0872 

rad 


-0.349 

0.489 

rad/sec 

“a 

-0.0872 

0.0872 

rad/sec 


This problem was then solved utilizing the proeedures deseribed in Section C of 
this chapter. The problem produced a feasible result and was then scaled by utilizing the 
following designer units to refine the solution; 

SCALES.X= 100; 

SCALES.Y= 100; 

SCALES.Z= 100; 

SCALES.V = 50; 

(5.15) 

SCALES.SIGMA= 1; 

SCALES.GAMMA = 0.1; 

SCALES.ALPHA= 1; 

SCALES.MU= 1; 
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b. Verification and Validation 

As with the maximum range OCP, the results for the minimum time problem were 
verified against the necessary conditions (omitted for brevity) and propagated through the 
dynamics equations for the problem, Equation (2.7), using the ODE45 solver in 
MATEAB. The results of the propagation were plotted against the DIDO solution for all 
of the states. The data are shown in Figures 66-69. 



Figure 66. Verification of Missile Trajectory for Minimum Time Maneuver 
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Figure 67. Verification of Missile Trajectory for Minimum Time Maneuver 



Figure 68. Verification of Velocity Vector for Minimum Time Maneuver 
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Figure 69. Verification of Flight Angles for Minimum Time Maneuver 

The propagated states for the minimum time problem shown in the Figures above 
were obtained by propagating the controls for the maneuver shown in Figure 70. 



Figure 70. Controls versus Time for Minimum Time Maneuver 
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c. Comparison of the Results 

Once the results are verified as in the previous seetion, they ean be eompared to 
the eonventional approaeh. When the trajeetories of the missile are plotted against eaeh 
other for both methods, shown in Figure 71, the differences between the two paths ean be 
seen. While the eonventional approaeh produeed a mostly linear response, the optimal 
eontrol solution eaused the missile to fly on a more parabolie trajeetory. 



Figure 71. Comparison of Missile Trajeetory for Conventional and Optimal 
Control Methods for Minimum Time Problem 


While the trajeetories show the two different paths that the missile takes to arrive 
at the same point, they do not provide any useful metrie for eomparing the two results. 
When the time of flight is plotted for both of the solutions. Figure 72, the effeets of the 
different trajeetories ean be seen. As mentioned earlier in this ehapter, the missile impaets 
the ground at a time of 32.3 seeonds when using the conventional proportional guidance 
based eontrol. When the missile is flown utilizing the optimal eontrol solution, the 
missile is able to reaeh the target in a time of 19.1 seeonds. The optimal eontrol solution 
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results in a reduction of about 40% of the flight time required for the conventional 
approach, which can provide an advantage for the warfighter. 



Figure 72. Time of Flight Comparison for Conventional and Optimal Control 

Methods for Minimum Time Problem 

The last comparison that was done for the minimum time problem was a 
comparison of the final velocity upon impact, shown in Figure 73. The conventional 
approach produces a final velocity of 130.2 m/s on impact, while the optimal control 
solution produced a final velocity of 401.0 m/s on impact. The optimal control solution 
velocity is greater than the conventional approach velocity by a factor of approximately 
3. Depending on the desired target, an increase in the maximum velocity may allow the 
missile to penetrate targets that had previously been impenetrable using conventional 
means. 

The terminal velocity obtained using optimal control is significantly larger than 
proportional navigation. To determine whether the increase is possible on the real missile 
(without violating other unmolded constraints), a more detailed simulation study must be 
performed, including a 6DOF simulation and analysis of structural loads. If the solution 
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developed using the reduced order model is feasible for flight implementation a 
significant tactical advantage may be possible simply by updating the missile flight paths. 



Figure 73. Comparison of Velocity versus Time for the 
Minimum Time Maneuver 


E. SUMMARY 

In this chapter optimal control theory was used to obtain a guidance solution to 
maximize the range of the fictitious missile under study. The bounds and conditions that 
were used for the problem setup as well as development of the necessary conditions for 
optimality were shown. The results of the solutions obtained using DIDO were presented 
and discussed. A minimum time maneuver of the fictitious missile was also obtained and 
the results of the solution obtained by DIDO were presented and discussed. The resulting 
maneuvers were compared to the conventional proportional navigation based solution and 
significant enhancements of the missile performance were observed using the optimal 
control-based approach. 
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VI. CONCLUSIONS AND FUTURE WORK 


A. CONCLUSIONS 

The intent of this thesis was to apply advanees in optimal eontrol theory, namely 
pseudospeetral optimal eontrol, to improve the range of a missile without ehanging any of 
the physieal eharaeteristies of the missile or its hardware. This software only approaeh 
eould be used to replaee eonventional guidanee logie, whieh has not ehanged in over 
50 years. Additionally, this thesis was intended to eontribute useful researeh that ean be 
built on for other applieations sueh as intereontinental ballistie missiles and related spaee 
systems. 

A redueed-order 3DOF model with simplifying assumptions was used to generate 
an optimal missile trajeetory that provided nearly 300% inerease in the maximum range 
of the missile over the eonventional proportional navigation-based approaeh. The newest 
variants of the medium-range air-to-ground missiles in the United States Order of Battle 
have a ranges that tend to fall from 80-100nm. An inerease in range in the terminal phase 
of flight by 8.75 nm is therefore a signifieant improvement in the reaeh of the missile. 

In a minimum time maneuver, it was shown that the time of flight eould be 
redueed by about 40% over eonventional proportional navigation. While the optimal 
eontrol solutions eould be obtained in a reasonable amount of time, the MATLAB based 
proeess eould not produee the solution in real time. Nonetheless, eomputing an optimal 
trajeetory prior to flight allows for the optimal trajeetory to be implemented by storing 
the appropriate eommands. Moreover, it may be possible to obtain real-time solutions if a 
speeialized eode is developed for these problems, whieh would inerease the utility of the 
optimal eontrol-based approaeh for missile guidanee. 

B. FUTURE WORK 

The work that was done in this thesis provides a good starting point for eontinued 

work on this problem. This thesis was limited to investigating only the mid-to-terminal 

phase of a missiles flight path, onee the roeket motor was burned out. Expanding the 

flight profile to inelude the other stages of flight will provide a more eomplete pieture on 
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the advantages that ean be obtained through use of optimal eontrol. The results that were 
found when examining just the final stage of flight, suggest that if the optimal eontrol is 
applied throughout the entire flight of the missile, the range of a missile eould be further 
inereased or other performanee benefits obtained. In addition, it would be helpful to 
explore the range of eonventional mid-eourse guidanee algorithms as eompared to 
optimal eontrol-based solutions. 

The problem examined in this thesis assumed the missile would intereept a 
stationary target. One of the obvious next steps is to apply optimal eontrol theory to a 
missile intereept problem involving a moving target. The prineiples of finding an optimal 
path to intereept a moving target ean not only impaet the guidanee of missiles but ean 
impaet the other area where proportional navigation type laws are used, sueh as aireraft 
intereepting moving targets or satellites performing elose proximity rendezvous. 
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